{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Low Mach approximation Navier-Stokes element formulation and automatic differentiation\n",
    "**Rubén Zorrilla**\n",
    "UPC BarcelonaTech, CIMNE\n",
    "\n",
    "## Governing equations\n",
    "This formulation comes from a simplification of the fully compressible Navier-Stokes equation and is suitable for the resolution of low speed strongly thermally coupled flows.\n",
    "Specifically, the formulation indroduces the low Mach limit simplification, which allows to split the total pressure into two parts, a thermodynamic part $p_{th}$ that is uniform in space and a mechanical (hydrodynamic) one $p$.\n",
    "It is important to note that $p$ is several orders of magnitude smaller than $p_{th}$ and is therefore omitted from the state and energy equations.\n",
    "In other words, the compressibility of the fluid comes from the thermal effects only.\n",
    "Hence, the problem to be considered is finding a velocity field $\\bm{u}$, mechanical pressure $p$, temperature field $T$ and thermodynamic pressure $p_{th}$ such that\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    \\frac{\\partial \\rho}{\\partial t} + \\nabla\\cdot\\left(\\rho\\bm{u}\\right) &= 0 \\\\\n",
    "    \\rho\\frac{\\partial\\bm{u}}{\\partial t} + \\rho\\bm{u}\\cdot\\nabla\\bm{u} - \\nabla\\cdot\\boldsymbol{\\tau} + \\nabla p + \\sigma(\\bm{u} - \\bm{u}_{s}) &= \\rho\\bm{g} \\\\\n",
    "    \\rho c_{p}\\frac{\\partial T}{\\partial t} + \\rho c_{p}\\bm{u}\\cdot\\nabla T - \\nabla\\cdot\\left(\\kappa\\nabla T\\right) - \\alpha T\\frac{d p_{th}}{dt} &= Q\n",
    "\\end{aligned}\n",
    "\\; ,\n",
    "\\end{equation*}\n",
    "being $c_{p}$ the specific heat coefficient at constant pressure, $\\kappa$ the thermal conductivity and $\\alpha$ the thermal expansion coefficient.\n",
    "The system must be closed by a state equation relating $\\rho$, $T$ and $p_{th}$.\n",
    "Assuming an ideal gas behavior, the state equation is $\\rho=p_{th}/RT$, with $R$ being the specific gas constant that relates heat capacity coefficients at constant pressure and constant volume as $R = c_{p} - c_{v}$.\n",
    "Introducing the heat capacity ratio $\\gamma = c_{p} / c_{v}$ allows to alternatively write the state equation for ideal gases as\n",
    "\\begin{equation*}\n",
    "    \\rho = \\frac{\\gamma p_{th}}{c_{p}(\\gamma - 1)T}\\; ,\n",
    "\\end{equation*}\n",
    "something that at the same time makes possible to obtain the following magnitudes as\n",
    "\\begin{equation*}\n",
    "    \\frac{\\partial \\rho}{\\partial t} = \\frac{\\gamma}{c_{p}(\\gamma - 1)}\\left(\\frac{1}{T}\\frac{d p_{th}}{dt}-\\frac{p_{th}}{T^{2}}\\frac{\\partial T}{\\partial t}\\right)\\; ,\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "    \\nabla \\rho = -\\frac{\\gamma p_{th}}{c_{p}(\\gamma - 1)T^{2}}\\nabla T\n",
    "\\end{equation*}\n",
    "and\n",
    "\\begin{equation*}\n",
    "    \\alpha = \\left.-\\frac{1}{\\rho}\\frac{\\partial \\rho}{\\partial T}\\right|_{p} = \\frac{1}{T}\\;.\n",
    "\\end{equation*}\n",
    "\n",
    "The determination of the thermodynamic pressure depends on the nature of the flow, which can be roughly divided in \\emph{open} and \\emph{closed} flows.\n",
    "For open flows (\\ie flows with Neumann boundary condition) $p_{th}$ is given at each time step by the boundary conditions.\n",
    "For closed flows (\\ie flows without Neumann boundary condition) of ideal gases with neither inflow nor outflow conditions $p_{th}$ can be obtained from the $T$ field as\n",
    "\\begin{equation*}\n",
    "    p_{th} = p_{{th}_{0}}\\frac{\\int_{\\Omega}\\frac{1}{T_{0}}}{\\int_{\\Omega}\\frac{1}{T}}\\; ,\n",
    "\\end{equation*}\n",
    "where the subscript $0$ indicates the initial magnitudes at $t=0$ and $\\Omega$ the problem domain.\n",
    "\n",
    "## Variational form\n",
    "The variational form (i.e. the functional to be automatically differenctiated) is obtained by multiplying previous equations by the corresponding pressure, velocity and temperature test functions $q$, $\\mathbf{v}$ and $w$. After integration by parts, this yields\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    - \\left(q, \\frac{\\partial \\rho}{\\partial t}\\right)_{\\Omega} &- \\left(q, \\nabla\\cdot(\\rho\\mathbf{u})\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\frac{\\partial\\mathbf{u}}{\\partial t}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}\\right)_{\\Omega} - \\left(\\nabla^{s}\\mathbf{v}, \\boldsymbol{\\tau}\\right)_{\\Omega} + \\left(\\nabla\\cdot\\mathbf{v}, p\\right)_{\\Omega} \\\\\n",
    "    &- \\left(\\mathbf{v}, \\sigma(\\mathbf{u} - {\\mathbf{u}_{s}})\\right)_{\\Omega} + \\langle\\mathbf{v}, \\mathbf{t}\\rangle_{\\Gamma} + \\left(\\mathbf{v}, \\rho\\mathbf{g}\\right)_{\\Omega} + \\left(w, Q\\right)_{\\Omega} - \\left(w, \\rho c_{p} \\frac{\\partial T}{\\partial t}\\right)_{\\Omega} \\\\\n",
    "    &- \\left(w, \\rho c_p \\mathbf{a}\\cdot\\nabla T\\right)_{\\Omega} - \\left(\\nabla w, \\kappa\\nabla T\\right)_{\\Omega} + \\langle w, \\kappa\\nabla T\\cdot \\mathbf{n}\\rangle_{\\Gamma} + \\left(w, \\alpha T \\frac{d p_{th}}{dt}\\right)_{\\Omega} = 0\\; ,\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "being $\\mathbf{a}$ the convective velocity, which computed from the linearised previous iteration velocity $\\mathbf{u}^{k-1}$ and the mesh velocity $\\mathbf{u}_{m}$ as\n",
    "\\begin{equation*}\n",
    "\\mathbf{a} = \\mathbf{u}^{k-1} - \\mathbf{u}_{m}\\; .\n",
    "\\end{equation*}\n",
    "\n",
    "## Variational MultiScales stabilization\n",
    "The scales separation reads\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    p &= p_{h} + \\tilde{p} \\\\\n",
    "    \\mathbf{u} &= \\mathbf{u}_{h} + \\tilde{\\mathbf{u}} \\\\\n",
    "    T &= T_{h} + \\tilde{T}\n",
    "\\end{aligned}\\; ,\n",
    "\\end{equation*}\n",
    "\n",
    "being $p_{h}$, $\\mathbf{u}_{h}$ and $T_{h}$ the finite element scales and $\\tilde{p}$, $\\tilde{\\mathbf{u}}$ and $\\tilde{T}$ the subgrid scales, which are modelled from the finite element residuals as\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "\\tilde{p} & := \\tau_{c} R_{c}(\\mathbf{u}_{h},T_{h}) \\\\\n",
    "\\tilde{\\mathbf{u}} & := \\tau_{m} \\mathbf{R}_{m}(p_{h}, \\mathbf{u}_{h}, T_{h}) \\\\\n",
    "\\tilde{T} & := \\tau_{e} R_{e}(\\mathbf{u}_{h}, T_{h}) \\\\\n",
    "\\end{aligned}\\; ,\n",
    "\\end{equation*}\n",
    "\n",
    "being $\\tau_{c}$, $\\tau_{m}$ and $\\tau_{e}$ the stabilization parameters which are equal to\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    \\tau_{c} &:= \\frac{\\mu}{\\rho} + \\frac{c_{2}}{c_{1}}h\\lVert u\\rVert +\\frac{c_{3}}{c_{1}}\\frac{\\sigma h}{\\rho} \\\\\n",
    "    \\tau_{u} &:= \\left(c_{1}\\frac{\\mu}{h^2} + c_{2}\\frac{\\rho\\lVert u\\rVert}{h} + c_{3} \\frac{\\sigma}{h}\\right)^{-1} \\\\\n",
    "    \\tau_{e} &:= \\left(c_{1}\\frac{\\kappa}{h^2} + c_{2}\\frac{\\rho c_{p}\\lVert u\\rVert}{h}\\right)^{-1}\n",
    "\\end{aligned} \\; .\n",
    "\\end{equation*}\n",
    "\n",
    "$c_{1}$, $c_{2}$, and $c_{3}$ are algorithmic constants defaulted in the header file to $4$, $2$, and $2$, respectively. $h$ the characteristic element size.\n",
    "\n",
    "The finite element residuals $R^{c}(\\mathbf{u}_{h}, T_{h})$, $\\mathbf{R}^{m}(p_{h}, \\mathbf{u}_{h}, T_{h})$ and $R^{e}(p_{h}, \\mathbf{u}_{h}, T_{h})$ are defined as\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    &R_{c}(\\mathbf{u}_{h}, T_{h}) = -\\frac{\\partial\\rho}{\\partial t} + \\alpha\\nabla T_{h}\\cdot\\mathbf{u}_{h} - \\rho\\nabla\\cdot\\mathbf{u}_{h}\\\\\n",
    "    &\\mathbf{R}_{m}(p_{h}, \\mathbf{u}_{h}, T_{h}) = \\rho\\mathbf{g} - \\rho\\frac{\\partial\\mathbf{u}_{h}}{\\partial t} - \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{h} - \\nabla p - \\sigma(\\mathbf{u}_{h} - \\mathbf{u}_{s})\\\\\n",
    "    &R_{e}(p_{h}, \\mathbf{u}_{h}, T_{h}) = Q - \\rho c_p \\frac{\\partial T}{\\partial t} - \\rho c_p \\mathbf{u}_{h}\\nabla T_{h} + \\alpha T_{h} \\frac{d p_{th}}{dt}\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "Once the subscales are defined, one can derive the stabilization terms to be added to the Galerkin functional. Hence, after doing the required integration by parts to remove the subscale derivatives and dropping the high order terms (the use of linear elements is assumed), one obtains the final functional to which the automatic differentiation is to be applied\n",
    "\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "    - \\left(q, \\frac{\\partial \\rho}{\\partial t}\\right)_{\\Omega} &- \\left(q, \\nabla\\rho\\cdot\\mathbf{u}_{h}\\right)_{\\Omega} - \\left(q, \\rho\\nabla\\cdot\\mathbf{u}_{h}\\right)_{\\Omega} + \\left(\\nabla q, \\rho\\tilde{\\mathbf{u}}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\frac{\\partial\\mathbf{u}_{h}}{\\partial t}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\rho\\mathbf{a}\\cdot\\nabla\\mathbf{u}_{h}\\right)_{\\Omega} \\\\\n",
    "    &+ \\left((\\nabla\\rho\\cdot\\mathbf{a})\\mathbf{v}, \\tilde{\\mathbf{u}}\\right)_{\\Omega} + \\left((\\rho\\nabla\\cdot\\mathbf{a})\\mathbf{v}, \\tilde{\\mathbf{u}} \\right)_{\\Omega} + \\left(\\rho\\mathbf{a}\\cdot\\nabla\\mathbf{v}, \\tilde{\\mathbf{u}}\\right)_{\\Omega} - \\left(\\nabla^{s}\\mathbf{v}, \\boldsymbol{\\tau}\\right)_{\\Omega} \\\\\n",
    "    &+ \\left(\\nabla\\cdot\\mathbf{v}, p_{h}\\right)_{\\Omega} + \\left(\\nabla\\cdot\\mathbf{v}, \\tilde{p}\\right)_{\\Omega} - \\left(\\mathbf{v}, \\sigma(\\mathbf{u}_{h} + \\tilde{\\mathbf{u}} - {\\mathbf{u}_{s}})\\right)_{\\Omega} + \\left(\\mathbf{v}, \\rho\\mathbf{g}\\right)_{\\Omega} \\\\\n",
    "    &+ \\left(w, Q\\right)_{\\Omega} - \\left(w, \\rho c_{p} \\frac{\\partial T_{h}}{\\partial t}\\right)_{\\Omega} - \\left(w, \\rho c_p \\mathbf{a}\\cdot\\nabla T_{h}\\right)_{\\Omega} + \\left((\\nabla w)\\rho c_p \\mathbf{a}, \\tilde{T}\\right)_{\\Omega} \\\\\n",
    "    &+ \\left(wc_p (\\nabla\\rho\\cdot\\mathbf{a}), \\tilde{T} \\right)_{\\Omega} + \\left(w\\rho c_p (\\nabla\\cdot\\mathbf{a}), \\tilde{T}\\right)_{\\Omega} - \\left(\\nabla w, \\kappa\\nabla T_h\\right)_{\\Omega} + \\left(w, \\alpha (T_h + \\tilde{T}) \\frac{d p_{th}}{dt}\\right)_{\\Omega} = 0\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "Note that the quasi-static subscales assumption is taken so all the subscales time derivatives are dropped. Also note that $\\rho$ (and its spatial and time derivatives) is expressed as a function of the unknown variables $T_{h}$ and $p_{th}$ by means of the equation of state above.\n",
    "\n",
    "## IMPLEMENTATION\n",
    "The remaining implementation is similar to that of the weakly compressible Navier-Stokes element.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
